******************************************************************
***   Programmer: C. Deal                             
***   Date Created: March 6, 2024
***   Replication file for Deal & Gonzales 2023 (Pediatrics)
******************************************************************

*Tell Stata where your data are located and any files should be saved

*Change to directory of replication folder:
global path "C:\Users\cdeal\Box\New Box Folder for RAs\Cameron Deal\Completed Projects\Sexual Minorities & Homelessness\replication"
cd "$path"

*Open your data for analysis
use "$path\2017_2019_yrbs_states", clear

***Primary Grouping Variables

*Homelessness Binary
destring qwheresleep, generate(homelessnum)
gen homelessbin=.
replace homelessbin = 0 if homelessnum==1
replace homelessbin = 1 if homelessnum>=2 & homelessnum<=7
label define homelessbin_label 0 "Not Homeless" 1 "Homeless"
label values homelessbin homelessbin_label
label var homelessbin "Homelessness"
tab homelessbin

*Sexuality Binary (LGBT+ or not)
destring q66, generate(sexualitynum)
gen sexualitybin = .
replace sexualitybin = 0 if sexualitynum==1  
replace sexualitybin = 1 if sexpart2==3 & sexualitynum==1
replace sexualitybin = 1 if sexualitynum>=2 & sexualitynum<=4
label define sexualitybin_label 0 "Heterosexual" 1 "LGBT+"
label values sexualitybin sexualitybin_label
label var sexualitybin "LGBT+ or Not"
tab sexualitybin

*Gender Minority Binary 
gen genderbin = .
replace genderbin = 1 if qntransgender==1 | qtransgender=="3"
replace genderbin = 0 if qtransgender=="1" 
label define genderbin 0 "Cisgender" 1 "Gender Minority"
label values genderbin genderbin
label var genderbin "Gender Minority Status"
tab genderbin

*Sex
gen male = 0
replace male = 1 if sex==2
replace male = 2 if sex==.
label define male_label 0 "Female" 1 "Male" 2 "Missing"
label values male male_label

*Sexual Orientation Categories
gen sexualitycat=.
replace sexualitycat = 1 if sexualitynum==1 & sexpart2==1
replace sexualitycat = 1 if sexualitynum==1 & sexpart2==2
replace sexualitycat = 2 if sexualitynum==2 
replace sexualitycat = 3 if sexualitynum==3
replace sexualitycat = 4 if sexualitynum==4  
replace sexualitycat = 5 if sexpart2==3 & sexualitynum==1
label define sexualitycat_label 1 "Heterosexual" 2 "Lesbian or Gay" 3 "Bisexual" 4 "Unsure" 5 "Heterosexual with Same-Sex Partners"
label values sexualitycat sexualitycat_label
label var sexualitycat "Sexual Orientation Categories"
tab sexualitycat

*Age Categories
gen agecat1 = .
replace agecat1 = 1 if age>=1 & age<=3
replace agecat1 = 2 if age>=4 & age<=5
replace agecat1 = 3 if age>=6 & age<=7
replace agecat1 = 4 if age==.
label define agecat1 1 "12 and under to 14" 2 "15 to 16" 3 "17 to 18 and older" 4 "Missing"
label values agecat1 agecat1
label var agecat1 "Age Categories"
tab agecat1

*Grade Level is already found in the original dataset as grade
replace grade = 5 if grade==.
label define grade_label 5 "Missing"
label values grade grade_label
tab grade

*Race is already found in the original dataset as race4
replace race4 = 5 if race4==.
label define race_label 5 "Missing"
label values race4 race_label
tab race4

*Homelessness Categories 1- for comparison to broader non-homeless population
gen homelesscat=.
replace homelesscat = 1 if homelessnum==1
replace homelesscat = 2 if homelessnum==2
replace homelesscat = 3 if homelessnum==3
replace homelesscat = 4 if homelessnum==4
replace homelesscat = 5 if homelessnum>=5 & homelessnum<=6
replace homelesscat = 6 if homelessnum==7
label define homelesscat_label 1 "Not Homeless" 2 "Non-Parental Home" 3 "Shelter" 4 "Hotel" 5 "Streets" 6 "Other"
label values homelesscat homelesscat_label
label var homelesscat "Homelessness Categories"
tab homelesscat

*Homelessness Categories 2- for comparison within non-homeless population
gen homelesscat2=.
replace homelesscat2 = 2 if homelessnum==2
replace homelesscat2 = 3 if homelessnum==3
replace homelesscat2 = 4 if homelessnum==4
replace homelesscat2 = 5 if homelessnum>=5 & homelessnum<=6
replace homelesscat2 = 6 if homelessnum==7
label define homelesscat2 2 "Non-Parental Home" 3 "Shelter" 4 "Hotel" 5 "Streets" 6 "Other"
label values homelesscat2 homelesscat2
label var homelesscat2 "Homelessness Categories"
tab homelesscat2

*Destring sitename
encode sitename, gen(sitename_factor)
tab sitename_factor

*Policy Friendliness- These scores were compiled by the Movement Advancement Project (https://www.lgbtmap.org/equality-maps) and reflect the overall friendliness of state policies toward LGBTQ populations
gen policy_score=.
replace policy_score = -0.5 if sitename== "Arkansas (AR)"
replace policy_score = 34.75 if sitename== "California (CA)"
replace policy_score = 33.5 if sitename== "Colorado (CO)"
replace policy_score = 34.5 if sitename== "Connecticut (CT)"
replace policy_score = 25.25 if sitename== "Delaware (DE)"
replace policy_score = 29.75 if sitename== "Hawaii (HI)"
replace policy_score = 30 if sitename== "Illinois (IL)"
replace policy_score = 7.75 if sitename== "Kentucky (KY)"
replace policy_score = 35 if sitename== "Maine (ME)"
replace policy_score = 26.75 if sitename== "Maryland (MD)"
replace policy_score = 13.5 if sitename== "Michigan (MI)"
replace policy_score = 25.25 if sitename== "New Hampshire (NH)"
replace policy_score = 26.5 if sitename== "New Mexico (NM)"
replace policy_score = 4.75 if sitename== "North Carolina (NC)"
replace policy_score = 7.75 if sitename== "North Dakota (ND)"
replace policy_score = 16 if sitename== "Pennsylvania (PA)"
replace policy_score = 32.5 if sitename== "Rhode Island (RI)"
replace policy_score = -0.5 if sitename== "South Carolina (SC)"
replace policy_score = 34 if sitename== "Vermont (VT)"
replace policy_score = 18.75 if sitename== "Virginia (VA)"
replace policy_score = 15.5 if sitename== "Wisconsin (WI)"
label var policy_score "LGBT-Friendly Policy Score"
tab policy_score

*Stigma1- Measured by percentage favoring LGBT acceptance- These values were taken from the Pew Research Center's Religious Landscape Study, which measured the percent of a state that said homosexuality should be accepted (https://www.pewresearch.org/religion/religious-landscape-study/)
gen stigma_percent1=.
replace stigma_percent1 = 45 if sitename== "Arkansas (AR)"
replace stigma_percent1 = 69 if sitename== "California (CA)"
replace stigma_percent1 = 67 if sitename== "Colorado (CO)"
replace stigma_percent1 = 76 if sitename== "Connecticut (CT)"
replace stigma_percent1 = 67 if sitename== "Delaware (DE)"
replace stigma_percent1 = 64 if sitename== "Hawaii (HI)"
replace stigma_percent1 = 64 if sitename== "Illinois (IL)"
replace stigma_percent1 = 44 if sitename== "Kentucky (KY)"
replace stigma_percent1 = 69 if sitename== "Maine (ME)"
replace stigma_percent1 = 66 if sitename== "Maryland (MD)"
replace stigma_percent1 = 62 if sitename== "Michigan (MI)"
replace stigma_percent1 = 71 if sitename== "New Hampshire (NH)"
replace stigma_percent1 = 59 if sitename== "New Mexico (NM)"
replace stigma_percent1 = 54 if sitename== "North Carolina (NC)"
replace stigma_percent1 = 57 if sitename== "North Dakota (ND)"
replace stigma_percent1 = 63 if sitename== "Pennsylvania (PA)"
replace stigma_percent1 = 72 if sitename== "Rhode Island (RI)"
replace stigma_percent1 = 51 if sitename== "South Carolina (SC)"
replace stigma_percent1 = 79 if sitename== "Vermont (VT)"
replace stigma_percent1 = 61 if sitename== "Virginia (VA)"
replace stigma_percent1 = 66 if sitename== "Wisconsin (WI)"
label var stigma_percent1 "LGBT Acceptance"
tab stigma_percent1

*Stigma2- Measured by percentage favoring LGBT acceptance- These values were taken from the Nationscape Study, which measured the percent of a state that said homosexuality should be accepted (author's calculations) (https://www.voterstudygroup.org/nationscape).
gen stigma_percent2=.
replace stigma_percent2 = .7899506 if sitename== "Arkansas (AR)"
replace stigma_percent2 = .8844646 if sitename== "California (CA)"
replace stigma_percent2 = .8887743 if sitename== "Colorado (CO)"
replace stigma_percent2 = .9161621  if sitename== "Connecticut (CT)"
replace stigma_percent2 = .8970887 if sitename== "Delaware (DE)"
replace stigma_percent2 = .8941629 if sitename== "Hawaii (HI)"
replace stigma_percent2 = .8859923 if sitename== "Illinois (IL)"
replace stigma_percent2 = .8206061 if sitename== "Kentucky (KY)"
replace stigma_percent2 = .8907871 if sitename== "Maine (ME)"
replace stigma_percent2 = .8943089 if sitename== "Maryland (MD)"
replace stigma_percent2 = .8807039  if sitename== "Michigan (MI)"
replace stigma_percent2 = .9306204 if sitename== "New Hampshire (NH)"
replace stigma_percent2 = .8764906 if sitename== "New Mexico (NM)"
replace stigma_percent2 = .8344134 if sitename== "North Carolina (NC)"
replace stigma_percent2 =  .8635236 if sitename== "North Dakota (ND)"
replace stigma_percent2 = .886658 if sitename== "Pennsylvania (PA)"
replace stigma_percent2 = .9259561 if sitename== "Rhode Island (RI)"
replace stigma_percent2 = .8418053 if sitename== "South Carolina (SC)"
replace stigma_percent2 =  .9246795 if sitename== "Vermont (VT)"
replace stigma_percent2 = .869445  if sitename== "Virginia (VA)"
replace stigma_percent2 = .8916968 if sitename== "Wisconsin (WI)"
replace stigma_percent2=100*stigma_percent2
label var stigma_percent2 "LGBT Acceptance (2)"
tab stigma_percent2



***Health Outcomes

**Mental Health

*Considered Suicide
gen considered_suicide=.
replace considered_suicide = 1 if q26=="1"
replace considered_suicide = 0 if q26=="2"
tab considered_suicide

*Attempted Suicide
gen attempted_suicide=.
replace attempted_suicide = 1 if qn28==1
replace attempted_suicide = 0 if qn28==2
tab attempted_suicide

*Planned Suicide
gen planned_suicide=.
replace planned_suicide = 1 if q27=="1"
replace planned_suicide = 0 if q27=="2"
tab planned_suicide

*Sad or Hopeless
gen sad_hopeless=.
replace sad_hopeless = 1 if q25=="1"
replace sad_hopeless = 0 if q25=="2"
tab sad_hopeless

*Suicide attempt that required medical treatment
gen suicide_treatment=.
replace suicide_treatment = 1 if qn29==1
replace suicide_treatment = 0 if qn29==2
tab suicide_treatment

*Mental Health Index
gen mental_index= sad_hopeless + planned_suicide + attempted_suicide + considered_suicide + suicide_treatment
tab mental_index

*Suicidal Activities Binary
gen suicide_bin=.
replace suicide_bin = 1 if considered_suicide==1 | attempted_suicide==1 | planned_suicide==1
replace suicide_bin = 0 if considered_suicide==0 & attempted_suicide==0 & planned_suicide==0
tab suicide_bin

**Aggression

*Physical Fight
gen fight=.
replace fight = 1 if qn17==1
replace fight = 0 if qn17==2
tab fight

*Carried a Weapon
gen weapon=.
replace weapon = 1 if qn12==1
replace weapon = 0 if qn12==2
tab weapon

*Aggression Index
gen aggression_index= fight + weapon
tab aggression_index

*Aggression Binary
gen aggression_bin=.
replace aggression_bin = 1 if aggression_index==1 | aggression_index==2
replace aggression_bin = 0 if aggression_index==0
tab aggression_bin

**Risky Sex Behaviors

*Sexually Active
gen sex_active=.
replace sex_active = 1 if qn61==1
replace sex_active = 0 if qn61==2
tab sex_active 

*Had sex with 4 or more partners
gen sex_4=.
replace sex_4 = 1 if qn60==1
replace sex_4 = 0 if qn60==2
tab sex_4

*Substances before sex
gen substance_sex=.
replace substance_sex = 1 if qn62==1
replace substance_sex = 0 if qn62==2
tab substance_sex

*Unprotected sex
gen unprotected_sex=.
replace unprotected_sex = 1 if q64=="2"
replace unprotected_sex = 0 if q64=="1" | q64=="3" | q64=="4" | q64=="5" | q64=="6" | q64=="7" | q64=="8" 
tab unprotected_sex

*Risky Sexual Behaviors Index
gen sex_index= sex_4 + sex_active + substance_sex + unprotected_sex
replace sex_index = 0 if sex_4!=1 & sex_active!=1 & substance_sex!=1 & unprotected_sex!=1
tab sex_index

*Risky Sexual Behaviors Binary
gen sex_bin=.
replace sex_bin = 1 if sex_4==1 | sex_active==1 | substance_sex==1 | unprotected_sex==1
replace sex_bin = 0 if sex_index==0
tab sex_bin

**Substance Use

*Cigarette smoking
gen cigarette=.
replace cigarette = 1 if q30=="1"
replace cigarette = 0 if q30=="2"
tab cigarette

*Alcohol Use
gen alcohol=.
replace alcohol = 1 if qn41==1
replace alcohol = 0 if qn41==2
tab alcohol

*Binge Drinking
gen binge_drinking=.
replace binge_drinking = 1 if qn42==1
replace binge_drinking = 0 if qn42==2
tab binge_drinking

*Marijuana
gen marijuana=.
replace marijuana = 1 if qn47==1
replace marijuana = 0 if qn47==2
tab marijuana

*Cocaine
gen cocaine=.
replace cocaine = 1 if qn50==1
replace cocaine = 0 if qn50==2
tab cocaine

*Smokeless Tobacco
gen no_smoke_tobacco=.
replace no_smoke_tobacco = 1 if qn37==1
replace no_smoke_tobacco = 0 if qn37==2
tab no_smoke_tobacco

*Cigarillos, Cigars, etc.
gen alt_smoking=.
replace alt_smoking = 1 if qn38==1
replace alt_smoking = 0 if qn38==2
tab alt_smoking

*Vaping
gen vape=.
replace vape = 1 if qn35==1
replace vape = 0 if qn35==2
tab vape

*Substance Use Index
gen substance_index= cigarette + alcohol + binge_drinking + marijuana + cocaine + no_smoke_tobacco + alt_smoking + vape
tab substance_index

*Tobacco Binary
gen tobacco_bin=.
replace tobacco_bin = 1 if qntb4==1 
replace tobacco_bin = 0 if qntb4==2
tab tobacco_bin

*Alcohol Binary- already coded as alcohol

*Other variables surrounding substance use (methamphetamines, heroin, unprescribed drugs) weren't asked in all of the sampled states and had small sample sizes



*****************************************
*** Identify Sample
*****************************************
*Adjust weights for pooled estimation (divide by the number of years)
gen one_year=0
replace one_year = 1 if sitename== "Colorado (CO)" | sitename== "Delaware (DE)" | sitename== "Illinois (IL)" | sitename== "Maine (ME)" | sitename== "Maryland (MD)" | sitename== "Michigan (MI)" | sitename== "North Dakota (ND)" | sitename== "Rhode Island (RI)" | sitename== "South Carolina (SC)" | sitename== "Vermont (VT)" | sitename== "Virginia (VA)" | sitename== "Wisconsin (WI)" 

gen weight_sample=.
replace weight_sample=(weight/2) if one_year==0
replace weight_sample=weight if one_year==1


*Tell Stata you are using survey weights
svyset psu [pweight=weight_sample], singleunit(certainty) strata(stratum)
set more off, perm

* Alternate method- best delineation between inclusion in study and non-inclusion
gen included1=0
replace included1 = 1 if homelessnum!=. & sexualitynum!=.
label define included1_label 0 "Not Included" 1 "Included"
label values included1 included1_label
label var included1 "Included in Study"
tab included1


***************************************************************************
***********   			     ANALYSIS    			     ******************
***************************************************************************
*Descriptive Statistics with `tabout' command
cd "$path"

*sample sizes:
tab sexualitybin if included1==1

*Weighted Estimates
svy, subpop(if included1==1 & sexualitybin==1): tab homelessbin, col
svy, subpop(if included1==1 & sexualitybin==0): tab homelessbin, col

*Proper Prevalence
svy: tabulate sexualitybin homelessbin, count format(%14.3gc)

***************************************************************************
***********   			     TABLE 1     			     ******************
***************************************************************************

*Actual Descriptive Statistics with tabout commands
tabout sexualitybin homelessbin using descriptives_lgbt.xls if included1==1, svy cell(row se) stats(chi2) nlab(count) f(4) sebnone replace
foreach var in agecat1 grade race4 sex homelesscat sexualitycat {
tabout `var' homelessbin using descriptives_lgbt.xls if included1==1 & sexualitybin==1, svy cell(col se) stats(chi2) nlab(count) f(4) sebnone append
}
tabout sexualitybin homelessbin using descriptives_straight.xls if included1==1, svy cell(row se) stats(chi2) nlab(count) f(4) sebnone replace
foreach var in agecat1 grade race4 sex homelesscat sexualitycat {
tabout `var' homelessbin using descriptives_straight.xls if included1==1 & sexualitybin==0, svy cell(col se) stats(chi2) nlab(count) f(4) sebnone append
}

***************************************************************************
***********   			     TABLE 3    			     ******************
***************************************************************************

*Health descriptive stats

foreach var in considered_suicide planned_suicide attempted_suicide suicide_treatment sad_hopeless aggression_bin sex_4 sex_active unprotected_sex substance_sex alcohol binge_drinking tobacco_bin marijuana cocaine {
tabout `var' homelessbin using healthdescriptives_lgbt.xls if included1==1 & sexualitybin==1, svy cell(col se) stats(chi2) nlab(count) f(4) sebnone append
}

foreach var in considered_suicide planned_suicide attempted_suicide suicide_treatment sad_hopeless aggression_bin sex_4 sex_active unprotected_sex substance_sex alcohol binge_drinking tobacco_bin marijuana cocaine {
tabout `var' homelessbin using healthdescriptives_straight.xls if included1==1 & sexualitybin==0, svy cell(col se) stats(chi2) nlab(count) f(4) sebnone append
}

*Assess Significance with regression
foreach var in considered_suicide planned_suicide attempted_suicide suicide_treatment sad_hopeless aggression_bin sex_4 sex_active unprotected_sex substance_sex alcohol binge_drinking tobacco_bin marijuana cocaine {
svy: logit `var' i.homelessbin if included1==1 & sexualitybin==1, or
}

foreach var in considered_suicide planned_suicide attempted_suicide suicide_treatment sad_hopeless aggression_bin sex_4 sex_active unprotected_sex substance_sex alcohol binge_drinking tobacco_bin marijuana cocaine {
svy: logit `var' i.sexualitybin if included1==1 & homelessbin==1, or
}


***************************************************************************
***********   			     FIGURE 1    			     ******************
***************************************************************************

global covariates "i.agecat1 i.race4 i.sex i.sitename_factor"
global covariates_nostates "i.agecat1 i.race4 i.sex"

*Tab with weighting- see how homeless categories interact with LGBT vs. heterosexual
svy: tab homelesscat sexualitybin if homelesscat!=1, col
*Chi2 test to see if relationship if difference between distributions is statistically significant
tab homelesscat sexualitybin if homelesscat!=1 [aw=weight_sample], col expected cchi2
display chi2tail(4, 61.4)
*Multinomial Logistic Regression
svy, subpop(if included1==1 & homelessbin==1): mlogit homelesscat2 i.sexualitybin $covariates, rrr


***************************************************************************
***********   			     TABLE 2    			     ******************
***************************************************************************
*Separate Logistic regression models for gay and straight subpopulations analyzing state policies vs. stigma
svy, subpop(if included1==1 & sexualitybin==0): logit homelessbin policy_score $covariates_nostates, or
svy, subpop(if included1==1 & sexualitybin==1): logit homelessbin policy_score $covariates_nostates, or
svy, subpop(if included1==1 & sexualitybin==0): logit homelessbin stigma_percent1 $covariates_nostates, or
svy, subpop(if included1==1 & sexualitybin==1): logit homelessbin stigma_percent1 $covariates_nostates, or
*Controlling for stigma percentage
svy, subpop(if included1==1 & sexualitybin==0): logit homelessbin policy_score $covariates_nostates stigma_percent1, or
svy, subpop(if included1==1 & sexualitybin==1): logit homelessbin policy_score $covariates_nostates stigma_percent1, or

*Marginal Effects for Sexual Minorities
svy, subpop(if included1==1 & sexualitybin==1): logit homelessbin policy_score stigma_percent1 $covariates_nostates, nolog 
margins, at(stigma_percent1 = (43.32212 63.32212)) atmeans post
test _b[1._at] = _b[2._at]
svy, subpop(if included1==1 & sexualitybin==1): logit homelessbin policy_score stigma_percent1 $covariates_nostates, nolog 
margins, at(policy_score = (12.82894 22.82894)) atmeans post
test _b[1._at] = _b[2._at]

*Marginal Effects for Heterosexual Youth
svy, subpop(if included1==1 & sexualitybin==0): logit homelessbin policy_score stigma_percent1 $covariates_nostates, nolog 
margins, at(stigma_percent1 = (43.32212 63.32212)) atmeans post
test _b[1._at] = _b[2._at]
svy, subpop(if included1==1 & sexualitybin==0): logit homelessbin policy_score stigma_percent1 $covariates_nostates, nolog 
margins, at(policy_score = (12.82894 22.82894)) atmeans post
test _b[1._at] = _b[2._at]

*Marginal Effects for Alternative stigma measure (Nationscape)
svy, subpop(if included1==1 & sexualitybin==1): logit homelessbin policy_score stigma_percent2 $covariates_nostates, nolog 
margins, at(stigma_percent2 = (78 98)) atmeans post
test _b[1._at] = _b[2._at]
svy, subpop(if included1==1 & sexualitybin==1): logit homelessbin policy_score stigma_percent2 $covariates_nostates, nolog 
margins, at(policy_score = (12.82894 22.82894)) atmeans post
test _b[1._at] = _b[2._at]

svy, subpop(if included1==1 & sexualitybin==0): logit homelessbin policy_score stigma_percent2 $covariates_nostates, nolog 
margins, at(stigma_percent2 = (78 98)) atmeans post
test _b[1._at] = _b[2._at]
svy, subpop(if included1==1 & sexualitybin==0): logit homelessbin policy_score stigma_percent2 $covariates_nostates, nolog 
margins, at(policy_score = (12.82894 22.82894)) atmeans post
test _b[1._at] = _b[2._at]


